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Abstract 

The results of a joint experimental and com- 
putational study on the flowfield over a periodically 
pitched NACA0012 airfoil, and the resultant lift varia- 
tion, are reported in this paper. The lift variation over a 
cycle of oscillation, and hence the lift hysteresis loop, is 
estimated from the velocity distribution in the wake 
measured or computed for successive phases of the cycle. 
Experimentally, the estimated lift hysteresis loops are 
compared with available data from the literature as well 
as with limited force balance measurements. Computa- 
tionally, the estimated lift variations are compared with 
the corresponding variation obtained from the surface 
pressure distribution. Four analytical formulations for the 
lift estimation from wake surveys are considered and 
relative successes of the four are discussed. 

1. Introduction 

As a part of a continuing research program on 
control of separated flows over airfoils and blades, 1,2 an 
experimental study of the phenomenon of dynamic stall 
over a periodically pitched airfoil was initiated about two 
years ago. The objective has been fundamental in scope: 
to advance the knowledge in the area, maintain in-house 
expertise, and aid in computational efforts. Initially, 
detailed phase-averaged flowfield measurements and flow 
visualization were carried out for specific cases of airfoil 
oscillation. These results have been summarized in Refs. 

3 and 4. 

During the analysis of the wake vorticity data, it 
occurred to us that the unsteady lift on the airfoil could 
be estimated from the vorticity flux shed into the wake. 


The analytical foundation and the various approximations 
for such estimation are deferred to a later section in the 
text. In short, the idea was based on estimating the 
change in circulation over time from the shed vorticity 
flux. Then by assuming a suitable convection velocity the 
change in the lift on the airfoil over the oscillation cycle 
could be estimated. The method produced lift hysteresis 
loops that had remarkable similarities with previous 
measurements. 5 

It was felt that the method deserved attention 
because the lift was obtained entirely from the wake 
survey. Determination of the forces on an oscillating 
airfoil is not an easy task. Force balance measurements 
can suffer from interference from structural resonances 
and static pressure distribution measurements can suffer 
from spatial resolution and sensor response limitations. 

Subsequently, the analytical foundation of the 
method was studied further. Alternate formulations, such 
as due to Theodorsen, 6,7 and in the format of the analysis 
of Wu, 8 were considered. The ’noncirculatory* component 
of the unsteady lift, due to the inertia of the fluid moving 
with the oscillating airfoil, was also considered following 
Theodorsen’s analysis. The details of these are discussed 
in section 2.3. Unfortunately, due to the small size of the 
airfoil and other experimental limitations, the lift hystere- 
sis loops could not be measured directly for comparison 
for cases involving dynamic stall. Only limited results 
could be obtained with a force balance for a case at a 
very low reduced frequency (k). The comparison of the 
balance measurement with the estimated lift hysteresis 
loop, however, was quite encouraging. 

It was felt during these deliberations that much 
insight could be gained from a computational study of the 
problem. One could calculate the unsteady lift from the 
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computed wake velocity field similarly as done in the 
experiment. Computationally, however, data on actual 
unsteady lift variation would be available from the cor- 
responding static pressure variation over the airfoil. Thus 
the validity of the estimated lift could be assessed by 
direct comparison. This led to a computational experi- 
ment, carried out by the third author. The computational 
results confirm the overall validity as well as potential 
deficiencies of the method under consideration. These 
results are described in the following. 

2 . Procedure and Analysis 

2. / Experimental procedure 

The experiments were carried out in a low speed 
wind tunnel the details of which have been described 
elsewhere. 2 The coordinate system with respect to the 
airfoil are schematically shown in Fig. 1. A two-dimen- 
sional model of a NACA0012 airfoil with chord, c = 

10.2 cm, and aspect ratio of 7.5, was mounted horizon- 
tally at the center (mid-height) of the test section. A 
pitching mechanism, described in detail in Ref. 3, was 
used to oscillate the airfoil about the one-quarter chord 
location. The oscillation amplitude and frequency could 
be adjusted continuously. An optical pick-up from the 
pitching mechanism was used to provide the reference 
signal for phase averaging. The experiments were con- 
ducted at chord Reynolds number, R c = 44,000. The 
angle of attack was varied as a - a tncim + a (l Sin(2ict/T), 
where T (= 1/0 is the period of oscillation. 

The flowfield measurements were carried out 
using a crossed hot-wire probe. The probe could be 
traversed in the streamwise (x) direction, and up and 
down in y, through an automated computer controlled 
traversing mechanism. All measurements reported are for 
the x-y plane at the mid-span location. The assumption of 
two-dimensionality is implicit in the investigation; data 
on the two-dimensionality of the flowfield have been 
presented in Ref. 3. 

As stated before, direct measurement of forces 
on the oscillating airfoil turned out to be difficult. Static 
pressure distribution measurement was not attempted due 
to a lack of availability of appropriate pressure trans- 
ducers small enough to be fitted in the airfoil model. A 
force balance, using load cells, 9 was used to measure 
steady lift variation with a. The same balance was tried 
in an effort to measure the unsteady forces for the oscil- 
lating case. The problem faced was harmonic distortion; 
typically a harmonic near the structural resonance would 
become large especially at a higher oscillation frequency. 
For the dynamic stall case with = 15°, the force 
changes were large as the airfoil went in and out of stall 
and the harmonic distortion was severe even at very low 
values of f. With = 0° and smaller amplitude (a, 

— 7.2°), the distortion was deemed minimal below an 
oscillation frequency of about 1 Hz. For this case the lift 


variation with a was measured and compared with the 
wake survey results as to be discussed in section 3.2. 

2. 2 Computational procedure 

The computational procedure is similar to that in 
an earlier study conducted by the third author and re- 
ported in Ref. 2. An upwind-biased, implicit, approxi- 
mate factorization algorithm which solves the thin-layer 
approximation to the two-dimensional Navier-Stokes 
equations is employed. 10 The algorithm is first -order ac- 
curate in time and second-order accurate in space. The 
computations are performed using a C-mesh. A quasi- 
one-dimensional characteristic analysis is used to expli- 
citly determine the far-field boundary conditions on the 
C-part of the mesh, while extrapolation is employed at 
the downstream boundary. On the airfoil, no-slip, adia- 
batic conditions along with zero normal pressure gradient 
are applied. For all results presented here, the Baldwin- 
Lomax algebraic turbulence model is used with transition 
location assumed to be at the leading edge. 

The unsteady motion of the airfoil is computa- 
tionally simulated by oscillating the grid as a rigid body 
about the one-quarter chord. 11 The converged steady state 
solution of the flowfield at is used as the starting 
point for the unsteady calculations. The computations are 
performed time accurately. For estimating the unsteady 
lift from the wake velocity field, the velocity and vor- 
ticity data are interpolated for stationary points in the 
flow field at a given downstream location in the wake. 

2.3 Analysis 

The analytical formulations are briefly described 
here. For further details the reader may consult Ref. 4, in 
which the experimental results can also be found in de- 
tail. The unsteady lift, L(t) for an oscillating airfoil can 
be divided into two components; ’non-circulatory’ L^t) 
and ’circulatory’ Leri). 7,12 The former component is due 
to the inertia of the fluid moved by the airfoil. This is 
dependent on k and can be negligible at small values of 
k. The ’circulatory’ component is due to the vortical flow 
arising from the airfoil surface. 

Non-circulatory part: Theodorsen provided an analysis 
for this component of the lift for pitching as well as 
plunging motion of a flat plate. 7 For pitching motion 
about one-quarter chord, as in the present investigation, 
the expression for the lift co-efficient can be written as 

k 2 

Q nc ( 0 - * ct # (k Cos2rcft — — Sin2irft). 

One observes that Cf NC increases as the square of k, and 
thus it becomes the dominant component at high values of 
k. It is linearly dependent on the amplitude (aj but 
independent of the mean angle (a^J. Finally, Cf NC at a 
given a is different between the upstroke and the down- 
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stroke, which yields a hystersis loop in its variation with 
a. 


Circulatory part: The subject of this paper is the circula- 
tory component which can be estimated from the vorticity 
flux in the wake. If one considers an impulsively started 
flow over a fixed airfoil, a ’starting vortex complex’ is 
created which convects away from the airfoil. Once the 
steady state is reached, the net amount of vorticity shed 
into the wake over a finite time is zero, and there is a 
constant circulation around the airfoil. The latter circula- 
tion for the ’bound vortex’, according to Kelvin’s theo- 
rem, is equal and opposite to that of the ’starting vortex 
complex’. The force acting on a pair of counter-rotating 
vortices of circulation +T and -T separated by a distance 
x is given by, 7 


Force = — (pTx) 
dt 

For the steady airfoil, the ’starting vortex complex’ 
moves away from the airfoil at the freestream velocity, so 
that dx/dt = U„, and this leads to the familiar Kutta- 
Joukowski theorem for steady lift, L = pUJ\ 

For the unsteady case of an oscillating airfoil, 
there is also a starting vortex complex, which should 
have a constant circulation. The circulation of the bound 
vortex, however, varies periodically with the oscillation. 
Within a finite time $t, the net amount of vorticity shed 
into the wake is nonzero, and again by Kelvin’s theorem, 
this amount should be equal and opposite to the change in 
the circulation of the bound vortex fiT occurring within 
the same time 6t. At any instant, the shed vorticity in the 
wake with circulation -5T and the corresponding change 
in the bound vortex t>T can be thought of forming as a 
counter-rotating vortex pair which are moving away from 
each other at a convection velocity U c ; U c is usually 
smaller than U*. Therefore, for the flow under con- 
sideration, the change in the lift in time it can be esti- 
mated as 


6L = pUST 

The change in the circulation iT can be found 
by considering the fixed path ABCD shown in Fig. I, 

For a sufficiently large path, it is reasonable to assume 
that all the vortical fluid is convected across the boundary 
CD only. For the two-dimensional, incompressible flow 
under consideration the time rate of change of circulation 
around the path ABCD is obtained, for example, from 
Eqn. 5.25 of Ref. 13 as 



dL 

dt 



Substituting £i(t) = / CD uo> z dy and integrating 
from time t=0 one obtains 


t 

L e (t) = - P U C / C,(t)dt + L c (0) (1) 

0 

The value at the beginning of the integration 
L c (0), cannot be determined from the vorticity flux and is 
assumed to be zero. However, this is just an additive 
constant and the integration, carried through a period of 
oscillation, should provide the same shape of the lift 
hysteresis loop regardless of the starting point for the 
integration. 

A caveat in Eqn. 1 is in the original formulation 
&L = pU c ar. The assumption that the unsteady forces 
are due to the interaction of the shed vortex and the 
corresponding change in the bound vortex neglects, for 
example, the interaction of the former with the bound 
vortex itself as well as with the starting vortex complex. 
The effects due to the interactions of the shed vortex with 
the latter two vortex systems, however, would mostly 
cancel each other, and thus Eqn. 1 may be a reasonable 
approximation. But the validity and accuracy have re- 
mained unclear. Note also that the formulation is equiva- 
lent to the application of the Kutta-Joukowski theorem to 
find the differential lift from the differential circulation, 
albeit using the convection velocity U c instead of U„ in 
the theorem. Of course, the convection velocity is not a 
clearly known quantity. It is expected to vary from case 
to case, e.g. when a a and a tnc4m are changed. Thus, the 
choice of U c in Eqn. 1 involves an arbitrariness. 

Alternate analyses for the unsteady lift, which 
lead to Eqns. 2, 3 and 4, are now described. Let us em- 
phasize here that all analyses considered have simplifica- 
tions and a foolproof method is not in sight. In an at- 
tempt to assess the validity of each equation, Eqns. 1-4 
will be applied to a given set of uc«) z (y,t) data. The resul- 
ting lift hysteresis loops will then be compared with 
available data, and with the actual lift obtained from pres- 
sure distributions in the case of the computation. 

The second equation for L c (t) is based on the 
flutter analysis of Theodorsen. 6,7 After certain simplifica- 
tions the expression for the lift co-efficient can be written 
in the present notations as 


L c(‘)= -pU c f 


C 

x + — 

2 


- 10T 




- & (2) 
2 2 


where cd 2 (== dv/dx-du/dy) is the spanwise component + L c (0), 

of vorticity. Therefore, one can write 

where as before, is equal to the vorticity flux Juco z dy 
which is a function of time. As in Eqn. 1, the term L c (0) 
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is unknown but is merely an additive constant. In order 
to evaluate the integral, at each time step a spatial distri- 
bution of vorticity is constructed from the temporal co z 
(y,t) data using a constant convection velocity U c (i.e., 
replacing x = -U c t). The integration is carried through 
the partial ’wavelength’ from the trailing edge and then 
over ten additional complete ’wavelengths’. (Increasing 
the number of complete wavelengths beyond the first 
partial wavelength did not make a significant difference 
in the result.) Thus, one finds that the primary difference 
between Eqns. 1 and 2 is the weighting factor inside the 
integral of the latter. The weighting factor is greater than 
unity close to the airfoil and quickly decreases to unity 
for x>c/2. Physically, it signifies that the vortices closer 
to the airfoil make a larger contribution to the net lift 
than those which are farther away. One also notes that the 
weighting factor in Eqn. 2 can be expressed in a binomial 
series as 


<* + i> 


(*+ 4) 2 - ( 4) 2 


(1 


+ 


¥ 


C/2 2 

x + c/2 



x + c/2 


). 


Thus the leading term in Eqn. 2 becomes the same as E- 
qn. 1 and the difference between the two lies in the con- 
tribution from the higher order terms. 

Equation 3 is adopted after Wu. 8 Starting with 
the Navier-Stokes equation Wu derived the expression for 

4(0 as 


MO = P /x« t dxdy 

This states that the force in the y direction (lift) is equal 
to the rate of change of the x-moment of all the vorticity 
in the flowfield. Note that the above equation is a gener- 
alized formulation of the relation, force = d/dt(pxT), 
used to obtain Eqn. 1. 

Thus, if the vorticity distribution over the entire 
flowfield were known the forces could be calculated 
accurately via the above equation. However, it would be 
practically impossible to measure all the vorticity, espe- 
cially around the oscillating airfoil. Thus, the following 
approximations are needed. From the measured distri- 
bution of <o 2 (y,t) at a given x, a spatial distribution is 
constructed by invoking the Taylor hypothesis, as was 
done for Eqn. 2. The unsteady lift, acting at mid-chord, 
is then approximated as 

t 

L c (0= -P 4 / (x + |)4(0dt + 4(0). (3) 

A -10T 2 

The fourth equation follows from the fact that 


the circulation around the airfoil at any instant is equal to 
the line integration of the velocity around the path ABCD 
(Fig. 1). For a large boundary ABCD this can be ap- 
proximated as, T(t) = f CD <v>(t) dy. Therefore, as 
done for Eqn. 1, one can write, 

4(t) = -pU c f <v>(0dy (4) 

CD 

Note that if <v>(t) represents the total phase 
averaged transverse velocity (including the long time 
average), the instantaneous total lift is evaluated by Eqn. 
4. 

Method of calculation: In the experiment, phase averaged 
axial velocities < u > , < v > and the spanwise compo- 
nent of vorticity < o> z > are used to evaluate all of the 
above equations. In the computation, the instantaneous 
values of these quantities are used. The full expression 
for the periodic variation of the lift coefficient C£ c is 
approximated, say from Eqn. 1, as 

U T l * 

CI c (0 = -2— — J j <u>*<o 2 >’dy *dt* . 

C 0 CD 

The superscript, *, represents non-dimensional ized quanti- 
ties (lift is nondimensionalized by l£pU„ 2 c, <o> 2 > by 
U./c, y by c and t by T). From the actual discrete mea- 
surements of <G) Z >*- and <u>* ij the above equation at 
any time step n+1, l<:n^NT, is evaluated as follows 

y j n NY 

Cl c (n + 1)= -2 £ <„>•„«■>,>* Ay/ At t '. 

c i=l j=l J 


As stated before, C£ c (l) is assumed to be zero. 

As discussed in Ref. 4, the measurements in the 
immediate vicinity of the airfoil trailing edge are marked 
by hot wire errors due to large flow angularity and occa- 
sional flow reversal. Thus all o> 2 (y,t) measurements are 
carried out at a downstream location, typically at x^/c 
= 0.3. This, however, introduces a time lag between the 
instants of measurement and the corresponding ’events’ 
taking place over the airfoil. This time lag is estimated 
as, -x_/U c , and accounted for in the calculation of the 
lift variation. 

An interesting condition arising from the re- 
quirement of finite lift on the airfoil is that the total 
change of lift over one complete period of oscillation 
should be zero. Therefore, the above calculation requires 

NT NY 

£ E <u> ’u <o> i > ‘ iJ A y/ At r = o- 

i=l J-l u 

Usually, due to measurement or computational errors this 
condition is not satisfied. This leads to Cf c (I) # 
Cf c (NT+l), i.e., an unclosed hysteresis loop. For 
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brevity, the deviation of this sum from zero is distributed 
over the entire cycle and only the resulting closed loops 
are presented. The ’closing error’ expressed as, 

NT NY 

E E < u> ’ y < < * ) z>’ y Ay " j At ’ , 

% error = i-1 -^ xlOO, 

NT NY 

E E r <«>•„<«,>* JAy-jAf, 

i«l j-1 y 

is listed in table 1 for a few cases considered in this 
paper. The denominator in the above expression repre- 
sents the sum of the absolute values of all vorticity shed 
in a cycle. This is referred to as the ’absolute vorticity 
flux’ which, in the experiments, was found to be approxi- 
mately a constant for a given value of and indepen- 
dent of k. The reader is referred to Ref, 4 for further 
discussion of these aspects, 

3. Experimental results 

Figure 2 shows a sequence of smoke-wire flow 
visualization photographs for a case involving dynamic 
stall, at various phases of the cycle. Frames (a) to (f) 
show phases when the angle of attack (a) is increasing 
(upstroke) and frames (g) to (j) show phases when a is 
decreasing (downstroke). As a increases, a clockwise 
vortex forms on the airfoil surface (frame d). This is the 
"dynamic stall vortex" as referred to by previous resear- 
chers. With further increase in a , the DSV moves to- 
wards the trailing edge. When it reaches the trailing edge, 
a counter-clockwise vortex starts to form near the trailing 
edge (frame 0- The counter-clockwise vortex, referred to 
as the trailing edge vortex (TEV), grows quickly under- 
neath the DSV (frames f and g) and lifts the latter from 
the airfoil upper surface. The DSV and the TEV combine 
to form a structure whose cross section looks like a 
mushroom. The ’mushroom’ structure evolves, moves 
upward and increases in size as it convects downstream 
(Frames h, i and j). In frame (i), at about 2 V4 chords 
from the trailing edge, its transverse extent is already 
very large and measures about 3 chords. After the pas- 
sage of the ’mushroom’ structure, frames (j) and (a) indi- 
cate the passage of a few smaller vortices before the cycle 
is repeated. 

While the DSV has been discussed in many 
previous papers on dynamic stall, the TEV and the ’mus- 
hroom’ structure have remained relatively unnoticed. 

Such structures were observed by only a few. 14,15 The 
intense TEV and the enormous ’mushroom’ structure 
could be quite significant in blade vortex interaction and 
aerodynamic noise generation, especially in configura- 
tions involving rows of blades. A computational study, 
using a multiple-scale turbulence model, was carried out 
recently for the conditions of Fig. 2 involving the dynam- 
ic stall. 16 Despite some differences with the experiment, 


the computation also yielded a similar sequence of events 
involving the DSV and the TEV. 

Detailed phase averaged flow field measurements 
were carried out for the flow condition of Fig, 2, and 
reported in Ref. 3. An example of the sets of data that 
led to the unsteady lift estimation is shown in Fig. 3(a). 
The temporal distribution of the phase averaged azimuthal 
component of vorticity, measured at x/c = 0.3, and over 
a complete cycle is shown in this figure. The axial and 
transverse velocity components, ensemble averaged typi- 
cally over 80 cycles, were measured at three closely 
spaced x-stations with x/c = 0.3 being the middle one. 
Central differencing provided 3<v>/3x(y,t) while the 
term 3<u>/3y(y,t) was evaluated by least squares 
fitting of the <u>(y,t) profiles. These two terms pro- 
vided the spanwise vorticity which was non-dimension- 
alized as, <cd z >* = <G) z >c/U„. 

While the pictures in Fig. 2 were obtained at an 
earlier date for k = 0.2, the somewhat lower k (=0.16) 
in Fig. 3(a) was chosen for use in evaluating Eqns. 1-4 
so that the estimated unsteady lift could be compared with 
available data from the literature. However, the sequence 
of events in the flow at these two values of k are very 
similar. A scrutiny of the data of Fig. 3(a) identifies the 
DSV and the TEV as the concentrated lumps of positive 
and negative vorticity around t/T = 0.5(a = 25°). (The 
sequence of events marked I, II, III and IV will be de- 
scribed later.) The vortical structures appear compressed 
in the temporal distribution because a full wavelength is 
captured in Fig. 3a while only a fraction of the wave- 
length is represented in the pictures of Fig. 2. The suc- 
cessive lumps of vorticity in Fig. 3(a), on the right of the 
TEV, represent a few additional vortices shed after the 
passage of the DSV-TEV pair. 

Figure 3(b) shows similar <*) z (y,t) data for a ttKm 
= 0° and a a = 7.2°, a case for which force balance 
measurements could be performed so that the lift varia- 
tion estimated from the vorticity data could be compared 
directly. One observes that because of the small a meMi and 
a a the vorticity distribution in this case is devoid of large 
concentrations and shows only a mild undulation. 

3.1 Lift variation estimated for the dynamic stall case 

The circulatory component of the lift coefficient, 
corresponding to the data set of Fig. 3(a), is shown in 
Figs. 4(a)-(d) as obtained by Eqns. 1-4, respectively. The 
data are shown as a function of a for a complete cycle. 
One finds that the estimates from Eqns. 1-3 are by and 
large comparable and differences occur mainly where 
there are steep variations. For example, the magnitude of 
the large drop in the lift around 25° is predicted differ- 
ently by the three equations. Since Eqn. 3 involves dif- 
ferentiation, the resulting curve appears somewhat ’jagg- 
ed’. The extent of the lift hysteresis loop, however, 
appears to be significantly ’underpredicted’ by Eqn. 4. 
Unfortunately, the relative accuracy of the four predic- 
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tions could not be judged directly as force balance mea- 
surement was not possible for this case (section 2.1). 

The circulatory component of Ct from Fig. 4(a) 
is added to the corresponding non-circulatory component 
(section 2.3) and the sum is shown in Fig. 5(a). The 
noncirculatory component is shown as the superimposed 
dashed curve. Since it is derived for inviscid flow about a 
flat plate the non-circulatory component should be re- 
garded as an approximation. However, it is small and the 
total unsteady Ct of Fig. 5(a) can now be compared with 
data from the literature. (As discussed before, the data in 
Fig. 5(a) only show departure from the steady lift, which 
is an unknown constant. In this figure the Ct values are 
referenced to the value at a = 5° where it is assumed to 
be zero). The data of McAlister et al., 5 obtained by static 
pressure distribution measurement at about the same value 
of k, are shown in Fig. 5(b). The overall features of the 
Ct variation in Fig. 5(a) can be seen to be similar to the 
data set of Fig. 5(b). The slope of the upper branch 
(between I and II) and the small anti -clockwise loop 
around a = 25° (between III and IV) in Fig. 5(a) are in 
reasonable agreement with the data of Fig. 5(b). The 
main difference occurs in the lower branch of the loop. 
But some differences are expected as the lift hysteresis 
loop is known to be sensitive to other flow parameters 
besides k. The Reynolds number R c was quite different 
between the two experiments (4.4xl0 4 in the present case 
as opposed to 4.8X10 5 in Ref, 5). The undulations on the 
lower branch, however, have been observed in other ex- 
periments. 17 Before continuing to evaluate the unsteady 
lift estimation technique, let us briefly discuss the ob- 
served variations in the lift in Fig. 5(a). 

Lift hystersis loop vis-a-vis measured vorticity: Obtaining 
the lift hysteresis loop from the vorticity data provided a 
unique opportunity to relate its various features with the 
vortical structures observed through the vorticity maps 
and the flow visualization. The variations in the lower 
branch of the lift hysteresis loop, as discussed in the 
foregoing, are believed to be real and due to the passage 
of the successive vortices following the DSV. As dis- 
cussed earlier, for the case under consideration, nearly all 
the positive (clockwise) vorticity generated from the 
airfoil suction surface accumulates to form the DSV 
during the upstroke (between points I and II in Figs. 3a 
and 5a). This is reflected in the wake as a depletion of 
positive vorticity. But the negative vorticity generated 
from the pressure surface is shed in the wake as usual. 
Qualitatively, a large negative vorticity in the wake is 
equivalent to a ’starting vortex’ and a large positive 
vorticity is equivalent to a ’stopping vortex’. When the 
former is shed, circulation around the airfoil as well as 
lift increases, while shedding of the latter causes a drop 
in the lift. Thus, between points I and II the airfoil lift 
increases and between points II and III, when the DSV 
containing positive vorticity is shed, the lift drops. The 


rebounding of the lift near the highest angle of attack 
during the downstroke (III to IV) is due to the shedding 
of the TEV which contains a concentration of negative 
vorticity. The undulations in the lower branch of the 
hysteresis loop occur due to the passage of a few more, 
relatively weaker positive and negative vortices following 
the DSV and the TEV (IV to I). 

3.2 Lift variation for a = 0° + 7.2°Sin2Kft 

The temporal distribution of vorticity for this 
case, at k = 0.028 (f = 0.57 Hz) and R c = 44,000, was 
shown in Fig. 3(b). The corresponding unsteady lift 
variation was measured with a force balance and also 
estimated using the above mentioned calculation proce- 
dures. The very low value of k was chosen to minimize 
harmonic distortions of the load cell signal (section 2.1). 
The force balance data are presented First followed by a 
comparative evaluation of the calculations. 

Force balance data: The unsteady lift variation measured 
by the force balance is shown in Fig. 6. The steady state 
lift variation, also measured by the same force balance is 
shown by the dashed line. The latter shows a kink around 
a = 0°, which, as discussed in the following, is believed 
to be due to laminar separation at this low operating 
Reynolds number. Such departure from linear variation 
due to laminar separation has been observed by others 
(e.g., Ref. 18). 

The unsteady measurements show a hysteresis 
loop even at this low oscillation frequency. The variation 
in the upper and the lower branches of the loop bear 
similarities with the steady state lift variation. At first 
sight, the hystersis loop is unexpected, since the dynamic 
stall phenomenon should not appear when the airfoil is 
oscillated within its static stall limit. 19 However, it is 
believed that laminar separation is responsible for the 
hystersis loop in much the same way as for the kink in 
the steady lift variation. For the steady airfoil, the flow 
remains separated on both surfaces around a = 0° result- 
ing in near zero lift. 18,9 Only when the angle of attack is 
increased (or decreased) sufficiently, does the flow reat- 
tach on the upper (or lower) surface resulting in the 
increase (or decrease) in lift. For the case of oscillation, 
the extent of the laminar separation on a given surface of 
the airfoil presumably depends on the direction of mo- 
tion. In other words, the extent of the separation at a 
given value of a on a given surface of the airfoil during 
upstroke is different from that occurring during down- 
stroke. This apparently causes the observed hystersis loop 
in the Ct curve. 

Estimated lift variation : Figures 7(a)-(d) show the lift 
hysteresis loops constructed from the data of Fig. 3(b) 
using Eqns. 1-4, respectively. The solid line represents 
the calculated circulatory part and in each figure this is 
plotted such that the mean Ct at a = 0° matches the 
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corresponding steady state value of the C£. The steady 
state lift variation measured by the force balance is shown 
as the dashed line. The non-circulatory component is 
negligible at this low value of k. Again, the lift hysteresis 
loops obtained by Eqns. 1-3 are found to be essentially 
similar and differences occur in the details. These loops 
are also very similar to the actual Cl variation of Fig. 6 
and the maximum and minimum amplitudes are well 
represented. The lift curve obtained by Eqn. 4, however, 
is clearly different. It is not completely clear but it ap- 
pears that a limited transverse distance in the integration 
is responsible for this underprediction. The computational 
results in the next section seem to confirm this. It should 
be emphasized here that the data of Fig. 7 were very 
sensitive to small changes, especially in oscillation fre- 
quency and in hot-wire calibration. The vorticity flux was 
small and thus accurate measurement was difficult; every 
time these data were retaken there was some difference in 
the lift hystersis loop. 

The experimental results for the small amplitude 
low k case provided the confidence that the unsteady lift 
measurement technique under consideration can yield 
results that are reasonably in agreement with the actual 
lift variation. Computational results presented in the next 
section add further credence to the technique. 

4 , Computational results 

The computations were carried out for the NACA- 
0012 airfoil at R c =44,000. In order to obtain more cycles 
of oscillation for a given number of iterations and other 
limitations a value of k- 0.3 and a freestream Mach 
number (M) of 0.3 were chosen for these computations. 
Furthermore, in order to avoid separated flows, for which 
the application of the Baldwin-Lomax turbulence model is 
questionable, a small amplitude oscillation condition of a 
= 2° + 2 0 sin2icft was chosen. The nonzero was 
chosen so that numerical errors would show up in the 
cumulative vorticity flux used to compute the unsteady 
lift. (A symmetric condition with a mean =0 produced 
cancellation of errors with a resultant ’closing error’ 
equal to zero). Note that at the chosen R c , laminar separa- 
tion would be expected in the experiment even at these 
low values of and a a . But this is not the case in the 
computation as boundary layer transition is assumed from 
the leading edge. As stated before, the idea here was to 
compute the lift variations using Eqns. 1-4 and compare 
those with the exact lift variation obtained from the static 
pressure (Cp) distribution, so that the validity of the four 
equations could be assessed. 

First, the adequacy of computational mesh density 
and domain extent were tested. Figure 8(a) shows the lift 
co-efficient variation with time, obtained from the (re- 
distribution, for three mesh densities. Time in this and 
the following figures is nondimensionalized by the chord 
and the speed of sound, and referenced to the instant 


when the unsteady calculations are initiated. These results 
show little variation, particularly on the two finer grids, 
indicating grid-convergent results for the lift history. 
Thus, the middle mesh density of 193x73 was deemed 
sufficient for the computations. Similarly, Fig. 8(b) 
demonstrates little difference in the results for computa- 
tional domains extending more than 8 chords from the 
airfoil. Most of the later computations are performed with 
the domain having 14c extent. Although not shown the 
lift coefficient variation with time obtained by Eqns. 1-4 
also exhibits convergent results for the aforementioned 
grid density and domain. The effect of time step on the 
C p -determined lift history was also investigated. Three 
different time steps were tried, for a total of 349, 698 
and 1396 steps (iterations) per cycle. All results were 
nearly identical to those shown in Fig. 8. Hence, 698 
steps per cycle were deemed to be sufficient. 

The lift variation obtained from the Cp-distri- 
bution (Fig. 8) is shown as a function of a in Fig. 9 for 
a complete cycle. The corresponding non-circulatory 
component calculated from the equation given in section 
2.3 is shown by the dashed line in Fig. 9. The lift varia- 
tion goes through hysteresis even for these small values 
of a mcan and a a . The hysteresis is mainly due to the non- 
circulatory component occurring at the relatively high k. 

It is also apparent that the calculated non-circulatory 
component does not fully account for the observed extent 
of the hysteresis loop. The discrepancy is believed to be 
due to simplifying assumptions in the equation for the 
non-circulatory component. For example, the airfoil 
thickness and the effect of the boundary layer are not 
taken into account in the formulation. 

While the lift variation obtained from C p -distri- 
bution, e.g., in Fig. 8, is "exact", the predictions from 
Eqns. 1-4 involve a phase lag depending on the measure- 
ment station. If such a phase lag is not accounted for 
properly it would drastically alter the shape of the lift 
hysteresis loop. For simplicity, however, the rest of the 
computational results are presented as a function of time 
only. Fig. 9 serves to provide an idea how these time 
variations would convert into the hysteresis loops. 

The Cl variations obtained by Eqn. 1, with data 
from different x-stations, are compared in Fig. 10 with 
the actual Cl variation. Note that the results of Eqn. 1 
represent the unsteady components, and for easy com- 
parison these are plotted such that the average of each 
matches the average of the actual lift variation. The 
results of Eqn. 1 are shown without any correction for 
the phase lag. Recall that in constructing the hysteresis 
loops from the experimental data (section 3.1) such a 
correction was done. The data were referenced to the 
trailing edge of the airfoil assuming a convection velocity 
U c . As expected, Fig. 10 shows a progressive phase lag 
with increasing distance of the measurement station. The 
phase lag, of course, arises due to the fact that there is a 
finite time required for the ’events’ occurring over the 
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airfoil to convert down to the measurement station. 
Nevertheless, it is amply clear from Fig. 10 that the 
amplitudes of the lift variation are reproduced very well 
by Eqn. 1. 

The times at which the peak in the Cl occurred 
at different x (Fig. 10) are shown in Fig. 11. A straight 
line can be reasonably fitted through the data. This in- 
dicates that it is reasonable to assume a constant convec- 
tion velocity, in conjunction with the equations, over the 
x-range covered. The slope of the straight line, however, 
yields a value of U c approximately equal to the freestream 
velocity (U„). Note that for the data presented in section 
3, based on experimental estimates, U C =0.6U„ and 
U C =0.6U„ were used for the 15° ±10° and the 0° ±7.2° 
cases, respectively. As stated before, it is not unexpected 
that U c varies from case to case, and this represents one 
of the arbitrary aspects in the method under consider- 
ation. However, U c - U„ is clearly appropriate for the 
case chosen in the computation, and this is what has been 
used for the rest of the computational results. 

The extrapolation of the straight line in Fig. 1 1 
also shows that when accounting for the phase lag the 
data should be referenced to a location approximately 
0.3c from the leading edge (x/c = -0.7). (The horizontal 
line represents the time when the C p -determined lift has 
the peak.) In other words, only when the data are shifted 
by a time, fit = (x+0.7c)/U c , would the comparison 
with the actual Cl variation be appropriate. It is not clear 
if this reference location would be the same for the 
dynamic stall case. However, if that were the case, the 
hysteresis loops constructed from the experiment would 
appear slightly different as in those cases the phase lag 
was corrected by referencing the data to the trailing edge. 

Figure 12(a) shows the C£ -variation obtained 
with Eqn. 1 for three different sampling rates as indi- 
cated. The results indicate that 50 samples/cycle is ade- 
quate for resolving the lift variation. This rate was used 
for the experimental data and is also used for the rest of 
the computational data. 

Since a grid system moving with the airfoil was 
used for the computations, the data had to be interpolated 
for stationary points in the wake at a given downstream 
station in order to emulate the experimental procedure. 

Fig. 12(b) shows lift variations obtained with Eqn. 1 for 
the indicated number of integration points. A given num- 
ber of integration points are spaced in the y-direction at a 
given x. The distance between the points is variable; it is 
the smallest at y = 0 and increases geometrically with 
distance. Clearly, for the transverse extent chosen (14c), 

22 intgration points misses the vortical region in the wake 
and results in near zero lift amplitudes. As the spatial 
resolution is improved from 194 points (corresponding to 
a minimum spacing of 0,01c) to 300 points (minimum 
spacing of 0.001c), convergence in the result is achieved. 
The 300 point integration is performed for the rest of the 
results. 


It should be mentioned here that the exercises 
done in Figs. 9 through 12 were repeated with the other 
three equations and the same inferences were made. The 
predictions of Eqns. 1-4 from the same set of data, at x/c 
= 0.5, are compared in Fig. 13. The transverse extent 
over which the integration was performed was varied for 
each case. For Eqns. 1-3 changes of the transverse extent 
did not make a difference in the predicted lift variation. 
Clearly, Eqns. 1 and 2 did well in predicting the lift 
variation. As in the experimental results, Eqn. 3, due to 
the differentiation involved, yielded a somewhat Jagged’ 
lift variation. For small transverse extent in the integra- 
tion, Eqn, 4, as in the experimental results, underpre- 
dicted the amplitudes. (The transverse extent in the exper- 
iment was approximately 1.5c). However, when the in- 
tegration distance is sufficiently large (> 16c), clearly 
Eqn. 4 also does just as well as Eqns. 1 and 2 in the 
prediction. 

Thus, the computational results for the chosen 
flow confirm that the wake survey method may be a 
viable approach for determination of unsteady lift on an 
airfoil. 

5. Concluding remarks 

A method of estimating the unsteady component 
of lift from wake velocity surveys is considered in this 
paper. The analytical foundations are discussed and four 
alternate equations with different approximations are 
considered. It is found that the lift hysteresis loops esti- 
mated with most of these equations compare well with 
limited force balance data as well as with data from the 
literature. The method is a novel one and could be of 
interest in similar experiments in the future as the lift 
hysteresis loop is obtained strictly from wake surveys 
without direct force or static pressure distribution mea- 
surements. 

The computational experiment carried out for a 
specific low amplitude case of airfoil oscillation supports 
the inferences made from the experimental data. These 
results confirm that the amplitudes of lift variation ob- 
tained by most of the equations are well in agreement 
with the actual lift variation. The latter is obtained from 
the Cp-distribution over the airfoil. These results also 
reaffirm that the convection velocity to be used in the 
method can vary from case to case. Even though a change 
in the convection velocity represents a change in a mul- 
tiplicative constant for the result, its choice represents a 
main arbitrariness in the method. The computational 
results showed that in order to account for the phase lag 
in the estimated lift variation, the data obtained at a 
certain downstream station should be referenced to a 
location approximately 0.3 chords from the leading edge. 

Of the four considered, Eqns. 1 and 2 are found 
to be equally successful in predicting the unsteady lift. 
Since Eqn. 1 is relatively simple and has the similarity 


8 


with the theorem for calculating steady lift from the 
circulation, this would be our recommendation. With this 
equation the ’circulatory’ component of the lift is esti- 
mated as L c = density x U c x cumulative vorticity flux 
shed by the airfoil from the beginning of an oscillation 
period, where U c is an appropriate constant convection 
velocity. From the computational results, it appears that 
this equation is relatively insensitive to sampling rate and 
integration extent, but a sufficient number of integration 
points are necessary to obtain accurate results. Equation 4 
should also be considered attractive, as it requires the 
simplest of measurements involving only one component 
of velocity in the wake. However, the transverse extent 
over which this measurement needs to be done in order to 
obtain convergent results is large. Equation 4, thus, may 
not be suitable for most wind tunnel experiments. 
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Table 1. ’Closing error’ for a few cases. 


j Method 

a 

k 

%c!osing 

error 

Computation 
x/c = 0.5 

2°± 2° 

.3 

-0.71 

Computation 
x/c = 1.5 

2°± 2° 

.3 

-0.86 

Experiment 

15°± 10° 

.16 

4.04 

Experiment 

0°± 7.2° 

.028 

1.6 



Fig. 1 Schematic of airfoil, co-ordinate system, and 
control volume for calculation of unsteady circulation; 
"u" and "d" denote increasing a (upstroke) and de- 
creasing a (downstroke). 
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Fig. 2 Smoke-wire flow visualization photographs at 
different phases of oscillation cycle; k = 0.2, a = 
15°+ 10°sin2nft. Approximate a for pictures a to j 
are 5°u, 14°u, 20°u, 22°u, 24°u, 25°u, 25°d, 20°d, 
16°d, 12°d. 
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Fig. 3 Temporal distributions of <a),>c/U.; x/c — 
0.3; R c = 44,000. Contour levels start at 0.5 or -0.5 
and are at interval of 0.5 (solid lines for positive and 
dashed lines for negative vortidty). (a) a = 15.3° + 
9.7°sin(2icft-rc/2), k - 0.16 (b) a = 0° 4- 7.2°sin- 
2*ft, k = 0.028. 



Fig. 4 Circulatory component of Cf versus a estimat 
ed from the data of Fig. 3(a); a = 15.3° + 
9.7°sin(2itft-ic/2), k = 0.16. (a) to (d) obtained by 
using Eqns. 1 to 4, respectively. 
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Fig, 5 C t versus a, for oscillation at a » 15° + 
10°sin2icft. (a) Unsteady component obtained by add- 
ing results of Eqn. 1 (Fig. 4a) with ’ non-circula tory 9 
component (dashed curve); k = 0.16, R c = 44,000, M 
= .019. (b) Total Ct versus a; k = 0.153, R c = 4.8 x 
10 5 , M = 0.036. 



Fig. 6 Total C t versus a for the case of Fig. 3(b): a 
= 0° + 7,2°sin2rcft, k = 0.028. Dashed line for 
steady state variation of Ct with a. 



a, deg 


Fig. 7 Circulatory component of Ct versus ct estimat 
ed from the data of Fig. 3(b); a = 0° + 7.2°sin2nft, 
k = 0.028. Figs, (a) to (d) obtained by using Eqns. 1 
to 4, respectively. 







Fig. 10 Cf versus time obtained by Eqn. 1, from com- 
puted data at different x/c as indicated. 


Fig. 8 Computed total Cf versus time obtained from 
Cp-distribution for a = 2° + 2°sin2n:ft; k = 0.3, R c 
= 44,000, M = 0.3. Effects of: (a) grid density, (b) 
computational domain extent. 




x/c 


Fig. 11 Time when peak in Cf occurred at different 
x/c corresponding to the data of Fig. 10. 


Fig. 9 Cf versus a corresponding to the data of Fig. 
8; dashed curve for ’non-circulatory’ component. 
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Fig. 12 Ct versus time obtained by Eqn. 1 from the 
computed data at x/c = 0.5; (a) effect of sampling 
rate, (b) effect of spatial resolution. 


o from Cp data 


4c extent 




time 


Fig. 13 Ct versus time obtained from computed data 
at x/c = 0.5. For each case, integration is performed 
over four transverse extents as indicated. Figs, (a) to 
(d) obtained by using Eqns. 1 to 4, respectively. 
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